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Abstract 

A semi-discrete fourth-order finite volume method is introduced to solve multi-dimensional mag¬ 
netohydrodynamic (MHD) problems. This method employs a dimension by dimension approach 
followed by one dimensional fourth-order centrally weighted essentially non-oscillatory reconstruc¬ 
tion (CWENO) to compute higher order point values of the physical quantities. Density, momen¬ 
tum and energy are discretized as volume averages. However, magnetic fields are discretized as 
area averages on the grid cell interfaces. The conservative semi-discrete finite volume scheme is 
combined with the constrained transport (CT) technique to enforce the divergence-free constraint 
on the magnetic field. Various ID, 2D and 3D tests are performed to confirm fourth order accuracy, 
shock capturing nature and divergence-free property of the scheme. 
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I. INTRODUCTION 


Many astrophysical phenomena that exhibit discontinuous solutions can be approximately 
described by the ideal magnetohydrodynamic (MHD) equations. Various shock capturing 
numerical methods have been proposed to solve astrophysical MHD problems [MB]. 

In this paper we introduce a semi-discrete finite volume scheme based on the computation 
of higher order local Lax-Friedrichs (LLF) fluxes (see for example |19[ 20]) at the grid cell 
interfaces. The LLF flux is an approximation to the Harten-Lax-van Leer (HLL) flux [21, 22] 
as it only requires the maximum propagation speed of information at the discontinuity which 
exists at the grid cell interface. The LLF fluxes are first computed at the centers of grid 
cell interfaces using ID fourth order CWENO reconstruction [23] following a dimension 
by dimension approach. Later an averaging is performed to obtain higher-order averaged 
LLF fluxes at the grid cell interfaces m ■ The CWENO approach is the heart of the 
present scheme as it straightforwardly allows us to compute a higher-order point value of 
the conserved quantities at the center of each grid cell interface. 

In this work, the constrained transport (CT) method [T] is used to evolve the magnetic 
held. This technique automatically maintains the solenoidality of the magnetic held up to 
the machine precision. In this approach the magnetic helds are usually discretized as face 
averages and electric helds as edge averages in the numerical grid cells. It is not straight 
forward to merge the CT technique into conservative hyperbolic equations due to different 
discretization of the physical variables. So far, several methods have been suggested on how 
to combine these two consistently [SHIM. In this paper, we use a very simple ID CWENO 
reconstruction [23] to compute higher order edge averaged electric held components from 
higher order LLF electric huxes. These edge averaged electric held components are used 
to evolve the magnetic held in the CT. At the end, a fourth order Runge-Kutta method 
[25] is used to evolve the semi-discrete scheme in time. We have performed several multi¬ 
dimensional tests to confirm the fourth order accuracy and the divergence free property 
(V.b = 0) of our scheme. 

The manuscript is organized as follows. Section II introduces the ideal MHD equations, 
the finite volume semi-discrete scheme for conservation laws, computation of averaged huxes 
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at the grid cell interfaces and the CT technique. Section III has been devoted to the ID 
fourth order CWENO reconstruction of all the physical quantities. In Section IV, we present 
the convergence of errors for ID, 2D and 3D Alfven wave tests at large amplitude to show 
fourth order accuracy in the nonlinear regime. We also present the convergence of errors for 
the 2D Orszag-Tang vortex problem as long as the solution is smooth. Moreover, the ID 
Shock tube, the Orszag-Tang vortex, the 2D blast wave and 2D Rotor problems are solved to 
demonstrate the shock capturing nature and divergence free property of the scheme. Section 
V contains a brief summary of the work presented in this manuscript. 

II. IDEAL MHD EQUATIONS AND SOLVER 

The ideal MHD equations, in conservative form, can be expressed as follows: 


dtp + V.(pv) = 0, 

(1) 

dt(pv) + V. [pvv + (p+ ^|b| 2 )/ - ^~ bb ] = °> 

(2) 

d t e + V.[(e + p+ ^|b| 2 )v - i(v.b)b] = 0, 

(3) 

<9jb + V x E = 0, 

(4) 


with the constraint 

V.b = 0. (5) 

Here b is the magnetic field, p is the density, e is the total energy density, v is the 
velocity, E = —v x b is the electric field. The thermal pressure p is computed from the 
ideal gas equation of state, 


where 7 is 
written as, 


p= ( 7 - l)[e- ^p|v | 2 - ^|b| 2 ]. 


( 6 ) 


the ratio of specific heats. The set of Eqs. 0-0 in compact form can be 


<9U <9F dG <9H 
dt dx dy dz 


(7) 
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where, 


U = 


p \ 

PV X 

P Vy 

PV Z 
e 

b x 

by 

\b z J 


( 


F = 


pv x 

pv 2 x +p + b 2 /87 r - b 2 /An 

pV X Vy - b X byj 47T 

pv x v z - b x b z /An 

v x {e +p + b' 2 /8n) - b x (v.b)/8n) 
0 

(V X by Vyb X ) 

(v z b x v x b z ) 


\ 


( 


, G = 


/ 


H 


P v y 

pv x v y - b x b y j An 

pv 2 + p + b 2 /87T - by 2 /An 
pV 2 y + P 

pv y v z - byb z /A7T 

v y (e +p + b 2 /87r) - fe y (v.b)/87r) 
(V X by Vyb x ) 

0 

0 yb z - v z by') 


\ 


pv z 

pv x v z - b x b z /An 
pvyv z - byb-jAn 
pv 2 z +p 

v z (e + p + b 2 /87r) — 6 2 (v.b)/87r) 
( v z b x v x bz) 


\ 


i v yb z 


V z by ) 


V 0 

All the symbols have their usual meanings. 
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( 8 ) 
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A. Semi-discrete scheme for conservative quantities 

Since density, energy and momentum components are discretized as volume averages, 
Eqs. 0-0 are solved by a semi-discrete method (see for example 0 ) as follows: 


7 — ttx _ ~nx r>V _ n,y 

dUij’k ' r i+l/2,j,k ' r i—l/2,j,k y jj+l/2,fc ^ i,j-l/2,k 

dt Ax Ay 

Here u l ,j.k is the volume average which is dehned as, 


+ 


n 


z 

i,j,k+ 1/2 


-m 


i,j,k— 1/2 




0. 


(9) 


r * z i,j,k-\-l /2 rVi,j+l/ 2 ,k f x i+l/ 2 ,j,k 




AxAyAz 


lZ(x, y, z ) dx d y d z . 


( 10 ) 


z i,j,k — 1/2 ^yi,j—l/2,k ** x i — l/2,j,k 


A 















The function 7Z(x, y, z ) is the higher order centrally weighted essentially non-oscillatory 
polynomial. The functions J rx , Q y and 7~L Z are the higher order averaged local Lax-Friedrichs 
(LLF) fluxes along x-, ?/-and ^-directions respectively. 


B. Computation of Averaged LLF Fluxes 


This subsection is devoted to the computation of the fluxes along the ^-direction as an 
example, the generalization to other directions is straightforward. We first reconstruct the 
ID CWENO polynomial along the x-direction from the given volume averages. Later this 
polynomial is used to compute the non-oscillatory point values centered at the cell interfaces 
along the x-direction. These reconstructed point values carry an area average information in 
the yz- plane. Considering them as cell averages, the ID reconstruction is performed along 
the ^-direction to obtain a point value in the cell center which still has an average information 
along the ^-direction. Again considering the reconstructed values as cell averages we perform 
the same ID reconstruction along the ^-direction to get the point value in the cell center. 
This is how we obtain the point values of the conserved quantities at the center of the grid 
cell interfaces along the x-direction. We can then compute corresponding point value fluxes 
f x and hence point value local Lax-Friedrichs (LLF) fluxes F x as follows : 


( 11 ) 


The quantities u E and u w denote right and left reconstructed point values of the conserved 
quantities at the center of the grid cell interfaces. These point values are obtained by ID 
CWENO reconstruction reconstruction [23] following dimension by dimension approach. 
The quantities ( f x ) E , ( f x ) w are corresponding point value fluxes of the conserved quantities 
at (x i+ i/ 2 ,yj,Zk). Here a x is the local maximum speeds of propagation which is estimated 


as, 


i+i/ 2 j,k = max{(|u x | + c f ) E jjk , ([w*| + 0}, 


where Cf is the fast magneto-sonic speed, 


c f = 


i h 2 i 1 / 2 

l(M + 4) + (M + c ?>) 2 -4c^) 1/2 ) 


( 12 ) 


(13) 
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Here c s = ( 7 p/p) 1 ^ 2 is the sound speed and ca = (\h\ 2 /Anp) 1 / 2 is the Alfven speed [7j. 
After we have computed the higher order point value LLF fluxes F x , the higher order 
averaged LLF fluxes T x can be easily computed using simple averaging along the y and 
^-directions [24] as, 


~rx _ t?x 1_ ( rpx _ 9 t?x _i_ j^x \ 

*'i+l/2,j,k ”i-\-l/2J,k ' \ "i-\-l/2,j—l,k i-\-l/2,j,k ' ”i-\-l/2,j-\-l,k) 


?x _ 9 rpx 1 t?x \ 

i+l/2,j,k—l ^ uJ -i-\-l/2,j,k ' ” i-\-l/2,j,k-\-l) 


The method, described above, requires all physical quantities to be volume averages. 
However, magnetic held components b Xl b y and b z are discretized as area averages on the 
yz, xz and xy-planes respectively. We therefore compute volume averaged magnetic held 
components (B x , B yj B z ) from their area averages ( b x , b y , b z ) as follows, 


1. Lagrange interpolation is used to obtain a fourth order accurate area average value 
at the center of the cell. It is to be noted that interpolation is only allowed along 
the direction where the held remains smooth and the direction is determined by the 
solenoiclality property of the magnetic held mm- 


2. Now since each grid cell has three face average values, one at the center and other two 
on the opposite faces, we are able to apply Simpson’s 1/3 rule to compute a fourth 
order accurate volume average. 


Thus, fourth order accurate volume average approximations for the magnetic held compo¬ 
nents are expressed as, 


B x i,j,k „ 


32 (90 b xi+ H2,j, k + 60 b X i-i/2,j, k 20 b xi+3/2 ^ k + 3b xi+5/2J , k 5b xi . 3/2J , k 


+b X i_i/2j : k + b X i + i/2j,k 


(14) 


Byi,j,k „ 


22 ( T 606yjj_]y2,fe 206yjj_|_3/2 ; fc T ''^byi.j | 5 / 2 .A' bb yl J 3 / 2 .A' 


+ byij- 1 / 2 ,k + byij+ 1 / 2 ,k 


(15) 
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- 5 zi,j,k s* 


32 1/2 206^,j,fc+3/2 “f - 36^2,^‘,£+5/2 ^^zi,j,k— 3/2 


+fyzi.7.fc-l/2 + b_ 


zi,j,k— 1/2 ~r ^zi Jt 7 ,/c+l /2 


(16) 


For the computation of fluxes along the ^-direction, and are being treated like 
other volume averaged variables (density, momentum etc.). However, ^-component of the 
magnetic held is special as it does not require any reconstruction along the ^-direction, but 
it does require the reconstructions along the y and ^-directions as it also has an average 
information in the yz-plane. 

Computation of averaged huxes in other directions is trivial due to the dimension by 
dimension approach. 


C. Constrained transport 

Since magnetic fields are discretized as area averages, one may, therefore, introduce finite 
volume discretization on Eq.Q after applying Stoke’s theorem as follows [I] : 


d 

dt 


bxi— 




E Z i- 


1/2,.7+1/2,fc 


E Z i- 


1/2,/—1/2,fc 


Ay 


+ 


Eyi-l/2,j,k+l/2 


Eyi-l/2,j,k-l/2 


Az 


(17) 


d _ Ezi+l/2,j-l/2,k ~ Ezi-l/2,j- 1/2, k E X i,j-l/2,k+l/2 ~ -E'xj,/-l/2,fc-l/2 

Tt yi ’ j ~ i /2 ’ fc “ Air Az 


(18) 


d u 

— Ori n 


E, 


dt Jzi ' j ' k ~ 1/2 Ax ' Ay 

The variables b x , b y , b z are the area averaged magnetic held components and E X1 E y , E z 
are the edge averaged approximations of the electric held components. 

From Eqs.([T7|)-([l9]), one may easily show that d(V .b)i t j t k/dt = 0 if the differentiation is 
carried out by using central diherences along the cell edges. Thus, the evolved magnetic 
held remains solenoidal if so initially. 


yi+l/2,j,k-l/2 


Eyi-l/2,j,k-l/2 E xi j + i/ 2,fc—1/2 


E ■ ■ 

-‘—'cr.r.'i — 


xi,j—1/2,k—1/2 


(19) 


7 










The edge averaged electric field components are obtained from the higher order LLF 
electric fluxes. These fluxes exist because Eq.Q can also be written in divergence form as 
follows, 


d t b + V. 


V 


0 E z 
~E Z 0 


E y —E x 


-■ E j 

E x 

0 
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= 0 , 


Thus LLF electric flux vector along the ^-direction reads as, 


( 20 ) 


E ' 7 


*+5 ,j,k 


1 

2 



f 0 1 

w 

( 0 ^ 

E - 


-E z 

+ 

-E z 



V Ey ) 

i+l,j,k 

\ Ey ) 






b.) 

w 


E - 


by 

— 

by 



Iv 

i+l,j,k 


i,j,k- 


( 21 ) 


Note here that the x-cornponent of the electric flux E xx is zero because b x fj k is equal to 
b.ci+\, : j.. k - The vectors E IF , E K are the left and right reconstructed point values which can be 
computed by using the definition of the electric field 


E = —v x b. 


( 22 ) 


and a x is already defined in Eq.(12). Once we have computed point value LLF electric 


fluxes along the x-direction, we apply the same averaging formula as in Eq.(14) to obtain 
the face (area) averaged LLF electric fluxes as follows, 


cx _ -c\x _ (T? x _ I t?x \ 

^i-\-l/2,j,k *- J i-\-l/2j,k ' r\ a l,fc ^ J ^- J i-\-l/2,j,k ' ^~ J i-\-l/2,j-\-l,k) 


24 

1 


j_ (T?x _ r )T? x I t^x \ 

' 24 ^^ 1 + 1 / 2 ,j,k —1 ' *~ J i-\-l/ 2 ,j,k-\-l) • 


Now, since all the components of the electric fluxes S x carry an average information 
along the y and ^-directions, ID reconstruction along any of those directions will return an 
edge average value along the other direction. Thus, considering — £ xy as cell average, if we 
reconstruct the point value at the cell edges along the ^-direction, we get a contribution to 
the edge averaged electric field along the ^-direction E z . Similarly, reconstruction of £ xz at 



















the cell edges along the ^-direction contributes to the edge averaged electric held along the 
^/-direction E y . 

Similarly, one can compute averaged LLF electric huxes along y and ^-directions and there 
contributions to the edge averaged electric held components. At the end we take a mean of 
all the contributions to compute fourth order edge averaged electric held components. 

Finally, Eqs.Q and are evolved in time using a fourth order Runge-Kutta 

method [25] , once the L.H.S. of these equations are computed. 

The spatial accuracy of this scheme is now completely determined by the order of spatial 
reconstruction of the physical quantities. We use ID fourth order CWENO reconstruction 
ra for the same. 

Although, the fourth order CWENO reconstruction for ID hyperbolic problems is well 
described by Levy et al. [23], we provide it here for the sake of completeness. 


III. ID FOURTH ORDER CWENO RECONSTRUCTION 

In each cell Ij, one has to reconstruct a quadratic polynomial Rj(x) which is a convex 
combination of three quadratic polynomials Pj_ i(x), Pj(x) and Pj+\ (x) such that, 

i+i j+ 1 

Rj(x)= ^2 w k p k (x ), where w k = 1, w k > 0, V k e (j - 1, j, j + 1). (23) 

k=j —1 k=j —1 

Rj (x) is reconstructed so as to satisfy all the three constraints accuracy, conservation and 
non-oscillatory. The coefficients of the polynomial P k (x ) are obtained uniquely by requiring 
it to conserve the cell averages u k ~ i, u k and u k+ 1 , where k G (j — 1 ,j,j + 1). Thus, each 
polynomial, P k (x), can be written as, 


P k (x) = u k - - 2 u k + u k _ i) + Uk+1 2 ^ k 1 {% ~ x k) 


(^fc+i 2 u k T u k — i) 
2Ax 2 


(x - x k ) 2 , k — j — 1 ,j,j + 1. 


(24) 


The nonlinear weights w k are obtained as before, i.e. 


w k = 


a k 


Oij —1 + Olj + i 


, where a k = 


Ck 


(e + IS k )P 


, Wk e(j - 1 ,j,j + 1). (25) 
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Here e, p are chosen to be 10 6 and 2 respectively and Cj_i = cy+i = 1/6, c, = 2/3 [23] . 
JS'fc are the smoothness indicators which are dehned as, 


is t = Y. 


<‘ x j+l/2 


(Ax)' 21 ix, Vfc e (j - i,j,j + 1). 


(26) 


1=1 


Zj-1/2 


Once we have reconstructed all the the polynomials {Pj~ i, Pj , Pj+i), smoothness indica¬ 


tors can easily be computed using Eq.(26) and hence nonlinear weights using Eq.(25). These 
weights are finally used to reconstruct non-oscillatory the polynomial Rj(x ) which provides 
fourth order accurate point values at the cell interfaces in smooth regions. 


IV. RESULTS FROM THE SIMULATIONS 

The fourth order accuracy of the proposed method has been verified by computing the 
convergence of errors in ID, 2D and 3D circularly polarized Alfven wave (CPAW) tests at 
large amplitude (A = 0.9) so as to retain the nonlinear behavior. We further confirm the 
fourth order accuracy in another 2D test — Orszag-Tang vortex problem when the solution 
is still smooth but highly nonlinear. 

In order to show the shock capturing behavior of the method we present the results for 
ID Brio-Wu test, the evolution of the Orszag-Tang vortex problem, the 2D MHD blast wave 
and 2D Rotor problem. The V.b errors in all the test are at the round-off level and do not 
grow in time. 

A. Accuracy 

In this subsection, we present convergence of the truncation error for the ID, 2D and 3D 
circularly polarized Alfven wave (CPAW) tests and for 2D Orszag-Tang vortex problem to 
validate the fourth order accuracy of the scheme. 

The norm of the error for all the convergence studies is computed as follows, 

1 N 

SE = ng'RM - ^ 1 - (27) 

i =1 

where A/, E{ are the reference and the numerical solutions as a function of grid resolution 
and AG is the number of grid points. 
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After computing the norms of the errors, we obtain the experimental order of convergence 
(. EOC ) using the formula, 


EOC(j) = 


\log(SE(NG(j)))\ - | log(6E(NG(j - 1)))| 


(28) 


\log(NG(j))\ — \log(NG(j — 1))| 
here j runs over the indices of the column vectors in the TABLES shown in the later 
subsections. 

CPAW tests — The initial conditions for ID, 2D and 3D CPAW tests are selected as 
follows: 


p — 1, v x = 0, v y = Asin[27rx], v z = Acos[27nr], 
p — 0.1, b x — 1 , by — Asin[27rx],. b z — Acos[27tx], 

[o,y x [o,y x [o,z z ] = [o, i] x [o, i] x [o,i]. (29) 


P = 1 , hs =-Acos [27 r(y + z)\i v y =-^=sm[2n(y + z)}, 
v z = ^=sm[2n(y + z)] t p = 0.1, b x = Acos[2n(y + z)\, 
by = 1 + -/= sin[ 27 r (y + z)], b z = 1 - -7= sm[ 2 n(y + z)\, 

[o,y x [o,y x [o ,i z \ = [0,1] x [o, i] x [0,1]. 


(30) 


P = 1, v x = - y cos[27t(x/Z x + y/ly + z/l z )\ - -^= sin[27r(x/Z x + y/l y + z/l z )\, 

ly l x V 2 

A 


A 


Vy = — cos[27T (x/l x + y/ly + z/l z )\ - -—-= 

LX ly\2j 


sin[27T (x/l x + y/ly + z/l z )\, 


A 


v z = —r = sm[ 2 n(x/l x + y/l y + z/l z )\, p = 0.1, 
v 2 

bx = ~ — 7 = - t- cos[27r [xjl x + y/ly + z/Z z )] - -— 7 = sin[27r(x/7 x + y/l y + z/Z z )], 
6„ = ^ + — cos[27t(x/^ + y/ly + z/Z z )] - -— 7 = sin[27r(a;/Z x + y/l y + z/Z z )], 

[y\ 2 


lyV2 la 


b z = 2= + — sin[2vr(a;// x + y/l y + z/Z z )], [0, Z x ] x [0, Z y ] x [0, Z z ] = [0, 2/v 7 ^] x [0, 2] x [0,1(31) 

It is clear from the initial conditions that for the ID case, CPAWs are propagating along 
the x-direct ion, for the 2D case, CPAW are propagating in the yz -plane at an angle n/4 to 
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the rc-axis and for the 3D case, CPAW propagate along the diagonal plane at an angle 7 t/ 4 
to the z-axis. 7 and A are chosen to be 5/3 and 0.9 respectively for all tests. Boundary 
conditions are periodic and the numerical solutions are obtained after one wave period. 


TABLE |TJ TABLE [TT] and TABLE 111 contain the convergence of errors for the ID, 2D 


and 3D CPAW tests respectively. 


TABLE I: 


Convergence of errors for the ID CPAW test 

j 

NG 

5E 

EOC 

1 

16 

3.4341404679648566E-003 

- 

2 

32 

1.479398592867942IE-004 

4.54 

3 

64 

6.7699572887503336E-006 

4.45 

4 

128 

3.518670464735423IE-007 

4.26 

5 

256 

2.0612437922231232E-008 

4.09 


TABLE II: 


Convergence of errors for the 2D CPAW test 

j 

NG 

SE 

EOC 

1 

16 2 

6.7880543545184406E-004 

- 

2 

32 2 

3.1515227458324933E-005 

4.43 

3 

64 2 

1.6701819661858828E-006 

4.24 

4 

128 2 

9.504393562687020IE-008 

4.14 

5 

256 2 

5.6911207021282479E-009 

4.06 


TABLE III: 


Convergence of errors for the 3D CPAW test 

j 

NG 

5E 

EOC 

1 

16 3 

7.2081207206123044E-003 

- 

2 

32 s 

2.6369454325840709E-004 

4.77 

3 

64 3 

9.2593960801853198E-006 

4.83 

4 

128 3 

3.5135732619242231E-007 

4.71 

5 

256 3 

1.7011058257482867E-008 

4.36 
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Orszag-Tang Vortex test — The initial condition for the Orszag-Tang Vortex test [7j are 
expressed as, 


P = 7 2 , v x = - sin (y), v y = sin(x), v z = 0, 
p = 7, b x = — sin(y), by = sin( 2 a;), b z = 0 , 
[ 0 ,U x [0 ,ly] x [ 0 , L] = [0,27r] x [ 0 , 2 vr] x [0,27r]. 


The ratio of specihc heats 7 is chosen to be 5/3. TABLE IV contains the convergence of 
errors at t — 0.5, where the reference solution is obtained at higher resolution (1056 2 ). 


TABLE IV: 


Convergence of errors for the Orszag-Tang vortex problem at t = 0.5 

j 

NG 

5E 

EOC 

1 

32 2 

8.3441012034921325E-003 

- 

2 

96 2 

6.8376844996276276E-005 

4.37 

3 

160 2 

5.8164163709434336E-006 

4.82 

4 

288 2 

4.0090947783788522E-007 

4.55 


TABLES [T IV confirm the fourth order accuracy of the scheme. 


B. Shock capturing behavior 

Brio-Wu Shock Tube test — This test is performed to demonstrate the shock capturing 
nature of the proposed method. The initial conditions for Brio-Wu Shock Tube test |2j are : 


0, v x , v y , v z , p, b x , by, b z ) = (1, 0,0, 0,1, 0.75,1,0), if x < l x /2 
0, v x , v y , v g , p, b x , b y , b z ) = (0.125, 0, 0, 0, 0.1, 0.75, -1, 0), if x > l x f 2. 

Boundary conditions are open along the ^-direction. The system length (l x ), number of grid 
points (n x ) and 7 are chosen to be 2, 1024 and 5/3 respectively. Figures [l][4] contain the 
spatial profiles of p, p, v y and b y respectively at t = 0.25. 

Orszag-Tang Vortex test — The initial conditions for this test are already described in 
the previous subsection. Here we provide the space-time evolution of the density (p) and 
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pressure at a resolution 256 2 to demonstrate the shock capturing nature of our scheme. 
Figures [5]{6] show the contour plots of the density (p) and pressure respectively at t = 0.5 
when the solution is smooth, however, Figures [7j{8] contain the snap shots of the same at 
t = 3.0 when the system exhibits discontinuities. 

2D MHD Blast wave test — The initial conditions for the 2D MHD Blast wave test H3 
are chosen as follows: 


(p, v x , v y , v z , p) = (1, 0, 0,, 0,10), if \J{x- l x / 2) 2 + {y- l y / 2) 2 <0.1 
(p, v x , v y , v z , p ) = (1,0,0,0,0.1), if \J{x- l x / 2) 2 + {y- 4/ 2 ) 2 > 0.1, 

b x = 1/V2, b y = 1 /a/2 , b z = 0, 

where [ 0 , l x ] x [ 0 , l y \ = [ 0 , 1 ] x [ 0 , 1 ]. 


The ratio of specific heats 7 is chosen to be 5/3 and grid resolution as 256 2 . Figures 
[9 10 contain the snap shots of the density (p) and pressure at t — 0.15 which clearly exhibit 
strong shock capturing nature of the present scheme. 

2D Rotor problem — MHD Rotor problem is well described in ref. (H, 5J [ 8 :, : H] • The com¬ 
putational domain is chosen to be unit square [ 0 , l x \ x [ 0 , l y \ = [ 0 , 1 ] x [ 0 , 1 ] and the resolution 
is 256 2 . The ratio of specific heats 7 is 1.4. The density p = 10 and the velocity components 
v x = -vo(y - 0.5)/r 0 , v y = v 0 (x - 0.5)/r 0 for r < r 0 , where r = sj(x- l x / 2) 2 + (y - l y /2) 2 , 
r 0 = 0.1 and u 0 = 1- For r > ry, p = 1, v x — 0 and v y = 0, where r\ = 0.115. How¬ 
ever, for r 0 < r < r 1; p = 1 + 9/, v x = —fv 0 (y — 0.5)/r, v y = fv 0 (x — 0.5)/r, where 
f — (r 1 — r)/(ri — ro). The pressure p and the magnetic field components are uniform; i.e., 
p = 1, b x = 5 /\/ 47 r, by = 0. Figures 11 12 contain the snap shots of the density (p) and 


pressure at t = 0 . 2 . 

In all the tests V.b remains conserved and does not grow in time. 


V. CONCLUSION 

We have successfully developed a very simple fourth order finite volume scheme to study 
astrophysical MHD problems. Accuracy, shock capturing nature and divergence free prop¬ 
erty are confirmed through various multi-dimensional tests. The method, presented in this 
paper, opens a possibility for implementing even higher order schemes simply by using higher 
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order ID reconstruction along with higher order averaging of fluxes and higher order time 
integrator. 
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FIG. 5: Contour plot of the density at t — 0.5 for the Orszag-Tang Vortex test 
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FIG. 6: Contour plot of the pressure at t — 0.5 for the Orszag-Tang Vortex test 
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FIG. 7: Contour plot of the density at t — 3.0 for the Orszag-Tang Vortex test 
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FIG. 8: Contour plot of the pressure at t — 3.0 for the Orszag-Tang Vortex test 


20 

























FIG. 9: Contour plot of the density at t = 0.15 for the 2D MHD blast wave 
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FIG. 10: Contour plot of the pressure at t = 0.15 for the 2D MHD blast wave 
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FIG. 11: Contour plot of the density at t — 0.2 for the 2D Rotor problem 
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FIG. 12: Contour plot of the pressure at t = 0.2 for the 2D Rotor problem 
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